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Thermoelastic Properties of Olivine and Wadsleyite 
(Fe^Mgi-^SiO^ Their Relationship to the 410 km Seismic 
Discontinuity. 



We combine density functional theory (DFT) within the 
local density approximation (LDA), the quasiharmonic ap- 
proximation (QHA), and a model of vibrational density of 
states (VDoS) to calculate elastic moduli and sound ve- 
locities of a— and (5— (Fe^Mgi-^SKIU (olivine and wad- 
sleyite), the most abundant minerals of the Earth's up- 
per mantle (UM) and upper transition zone (TZ). Com- 
parison with experimental values at room-temperature and 
high pressure or ambient-pressure and high temperature 
show good agreement with our first-principles findings. 
Using our results, we investigate the discontinuities in 
elastic moduli and velocities associated with the a — > 
f3— (Fe I ,Mgi_ a; )2Si04 transformation at pressures and tem- 
peratures relevant to the 410 km seismic discontinuity. We 
find the compressional velocity contrast to be smaller than 
the shear velocity contrast, in agreement with the prelimi- 
nary reference earth model (PREM). 



1. Introduction 

Wadsleyite (/3— phase) is a high-pressure polymorph of 
olivine (a— phase), (Mgi_ x ,Fe a; )2Si04. These minerals are 
the main constituents of the Earth's upper mantle (UM) 
and upper transition zone (TZ) [Rinwood, 1975; Irifune and 
Ringwood, 1987; Putnis, 1992]. At ~13.5 GPa, the trans- 
formation from olivine to wadsleyite happens near 1600 K 
[Akaogi et al, 1989; Katsura and Ito, 1989]. This transfor- 
mation is associated with the 410 km discontinuity in seismic 
velocities in the Earth. The interpretation of seismic mod- 
els in terms of mineralogy and chemical composition requires 
knowledge of bulk (K) and shear (G) moduli as functions of 
pressure and temperature of constituent minerals. Elastic 
properties of iron-bearing olivine have been studied exper- 
imentally at ambient temperature using impulsively stimu- 
lated laser scattering (ISLS) to 17 GPa [Abramson et 
1997 
1998 



Brillouin spectroscopy (BS) to 32 GPa [Zha et al, 
ultrasonic interferometry (UI) at elevated pressure 
and temperature [see e.g. Liu et al, 2005], and resonant ul- 
trasound spectroscopy (RUS) to 1500 K [Isaak, 1992]. Elas- 
tic measurements on iron-bearing wadsleyite have been more 
limited in temperature and pressure range than for olivine. 
Data have been obtained using UI at high pressure and high 
temperature [Li and Liebermann, 2000; Liu et al, 2009], 
RUS and resonant sphere techniques at ambient pressure 
and high temperature [Mayama et al., 2004; Isaak et al., 
2010]. 

Computational predictions of iron-bearing a— and 
/3— phases have been reported only for static elasticity at 
high pressure [Nunez Valdez et al, 2010, 2011; Stackhouse 
et al, 2010]. High temperature first-principles calculations 
employing the quasiharmonic approximation (QHA), valid 



up to ~ 2/3 of the melting temperature, or molecular dy- 
namics (MD) methods, fundamental near and above melting 
temperatures, complement each other and have been used to 
obtain elastic moduli. These calculations are, however, re- 
stricted in number since they are computationally demand- 
ing [Karki et al, 1999; Wentzcovitch et al, 2004]. QHA cal- 
culations required vibrational density of states (VDoS) for 
strained atomic configurations at several pressures, that is, 
about 1000 parallel jobs [Da Silva et al, 2007; Da Silveira et 
al, 2008, 2011]. In this paper we use a new analytical and 
computational procedure [Wit and Wentzcovitch, 2011] to 
calculate K, G and compressional (Vp), shear (Vs), and bulk 
(V<s>) velocities of a — fi— (Mgi_ a; ,Fe :E )2Si04 phases. So far 
this method has only been tested on MgO and a— Mg2Si04 
at room pressure and high temperature. This approach uses 
only static elastic constants and phonon density of states for 
unstrained configurations, therefore reducing the amount of 
computational time, resources, and especially human effort, 
by one to two orders of magnitude. We then address con- 
trasts across the a — > (3 (Fe a; ,Mgi_ a ;)2Si04 transition near 
conditions of the 410 km seismic discontinuity. 

2. Theory 

The computational approach is based on density func- 
tional theory (DFT) [Hohenberg and Kohn, 1964; Kohn and 
Sham, 1965] within the local density approximation (LDA) 
[Ceperley and Alder, 1980]. At arbitrary pressures (P's), 
orthorhombic equilibrium structures of olivine (space group 
Pbnm, 28 atoms/unit cell) and wadsleyite (space group 
Imma, 28 atoms/primitive cell) were found using pseudopo- 
tentials, energy cutoffs, and k-point samplings as reported 
by Nunez Valdez et al. [2010, 2011]. Interatomic forces and 
stresses were smaller than 10 -4 Ry/a.u. Dynamical matrices 
were calculated using density functional perturbation theory 
(DFPT) [Baroni et al, 2001] via quantum-ESPRESSO [Gi- 
annozzi et al., 2009]. At each pressure, a dynamical matrix 
was obtained on a 2 x 2 x 2 q-point mesh for one atomic 
configuration only. In principle about 10 other configura- 
tions should be used as well, but here we are interested in 
elastic, not thermodynamic properties, and on the depen- 
dence of averages of phonon frequencies with strains, thus 
the current approximation is sufficiently accurate [Wu and 
Wentzcovitch, 2011; Nunez Valdez et al., 2010]. Force con- 
stants were extracted and interpolated on a 12 x 12 x 12 
regular q-point mesh to produce VDoS. 

Then we exploit the dependence of phonon frequencies 
on anisotropic composition to determine the thermal con- 
tribution to the Helmholtz free energy F, within the QHA 
[Wallace, 1972], that is, 

F (e, V, T) = U st (e, V) + ^ V) + 



+fc B T^ln |l - exp 



fi^q,m(e, V) 

k B T 



(1) 
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where q is the phonon wave vector, m is the normal mode 
index, T is temperature, U s t is the static internal energy 
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Figure 1. (Color online) Pressure and temperature dependence of elastic moduli K and G (calculated as Voigt-Reuss-Hill 
averages [Watt, 1979]), and velocities Vp, Vs, and V<s> for Fe-free olivine-wadsleyite (a, c, e, g) and Fe-bearing olivine- 
wadsleyite (b, d, f, h). First principles calculations within LDA (lines) are compared to available experimental data 
(symbols). Low-pressure-high temperature trends (dash lines) are outside the validity of the QHA. 



at equilibrium volume V under isotropic pressure P and in- 
finitesimal strain e, fi and kp are Planck and Boltzmann con- 
stants, respectively. Isothermal elastic constants are given 
by CT jkl = [^£1]^ with G = F + PV, the Gibbs 

energy and i, j, k, I = 1, . . . , 3. To convert to adiabatic elas- 
tic constants, we use Cf jkl = Cf jkl + ^SijSu, 

(Cv =heat capacity at constant V, and S ^entropy). Cjj ki 
can be expressed as functions of strain and mode Griinesein 
parameters. Assuming that the angular distribution of 
phonons is isotropic, isothermal elastic constants can be cal- 
culated without performing phonon calculations for strained 
configurations [Wu and Wentzcovitch, 2011]. This approx- 
imation is equivalent to assuming that thermal pressure is 
isotropic. Even though this assumption is not completely ac- 
curate [Carrier et al, 2007], the method is a good approx- 
imation within experimental uncertainties [Carrier et al, 
2008]. Static elastic constants previously computed [Nunez 
Valdez et al., 2010, 2011] are also used. 

3. Results and Discussion 

To the best of our knowledge, this is the first DFT- 
based report of aggregate properties of Fe-bearing olivine 
and wadsleyite at pressures and temperatures relevant to 
the UM and TZ. Fig. (1) shows the calculated pressure de- 
pendence of K, G, V P = ^j{K+^G)/p, V s = ^/G/p, and 

V<s> = \J K/p, with p ^density, of olivine and wadsleyite at 
several temperatures compared to experimental data. The 
general agreement between our predictions and measure- 
ments at room temperature and 1000K is excellent for x — 
[Zha et al., 1996; Isaak et al, 1989; Li et al., 1996; Zha et al, 
1997; Isaak et al, 2007] within the valid regime of the QHA 
[Yu et al., 2008]. For x = 0.125 our trends compare very well 
with experiments in the range x = 0.08 — 0.013 [Abramson 
et al., 1997; Zha et al., 1998; Isaak, 1992; Sinogeikm et al, 



1998; Li and Lieberrnann, 2000; Liu et al., 2009]. K and G 
increase with pressure, but the pressure dependence of K 
is stronger. Temperature decreases both K and G but the 
former is the most affected, this effect propagates to Vp and 
Vg. for both x = and x = 0.125, Fig. (1). 

Fig. (2) displays aggregate properties of a— and 
f3— phases as functions of temperature at ambient pressure. 
The prediction power of our DFT-based calculations is out- 
standing within the limit of the QHA for Mg-end mem- 
ber and Fe-bearing a— phases, [see Figs. (2a and c)], even 
though the iron content in experimental data is somewhat 
different from x = 0.125 [Isaak et al, 1989; Isaak, 1992]. At 
ambient conditions, our results for K and G agree quite well 
with experimental measurements for a— and /3— Mg2Si04 
phases [on-line, 2011]. We find that iron content increases 
K and decreases G for both phases. Experimentally this 
effect is small in the case of olivine, Figs. (2a) and (2b). 
For wadsleyite, measurements show that K is even more 
independent of Fe content within experimental uncertainty, 
while G seems to decrease, in agreement with our findings. 
Furthermore, looking at the temperature dependence of j3- 
K at ambient pressure [Fig. (2b)] we see that our predicted 



Table 1. Predicted contrasts in % across the a — > fi transi- 
tion at 13.5 GPa. 
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Figure 2. (Color online) Temperature dependence of K, G, Vp, Vs, and V<s> for Mg-end members olivine-wadsleyite (a, 
c) and Fe-bearing olivine-wadsleyite (b, d). First-principles calculations within LDA (lines) are compared to available 
experimental data (symbols) at P=0 GPa. 
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Figure 3. (Color online) Density, elastic, and velocity 
contrasts (lines) compared to laboratory and seismic data 
across the Fe-free (a) and Fe-bearing (b) a — > ft transi- 
tion. Y2008 Yu et al. [2008], SF Shearer and Flanagan 
[1999], PREM Dziewonski and Anderson [1981]. 



trend for x = 0.125 falls in between two experimental sets 
of data for x — 0.08 and x — 0.09, with measurements for 
x — 0.09 [Mayama et al, 2004] being smaller than those for 
x — 0.08 [Isaak et al., 2010]. These observations show that 
it is uncertain from experiments whether iron in the range 
x — 0.00 — 0.10 causes K to increase or decrease since its ef- 
fect on K is so small. In contrast, calculations show clearly 
the dependence of K on x. 

Acoustic velocities of a— and /3— phases decrease with 
increasing x at all pressures and temperatures, [see Fig. 



(lc,d,g, and h) and Fig. (2c and d). Our results for both 
phases with x — are in very good agreement with ex- 
perimental values, Fig. (lc and g). For x = 0.125, our 
velocities are smaller than measurements in samples with x 
in the range 0.08 — 0.12, Fig. (Id and h). Accurate data on 
elasticity of a— and /3— Mg2Si04 at UM and TZ conditions 
are critical for investigating the role of the a— to /3 trans- 
formation on the 410 km seismic discontinuity. Predicted 
contrasts, A, for a property M across id-t/J transforma- 
tion, AM = 2 {M x ,g — M x , a ) I (M x , a + M Xt p) x 100, indi- 
cate major changes in elastic moduli and velocities in these 
two phases near 410 km depth conditions, Table (1). Fig. 
(3) shows how contrasts slightly decrease with pressure and 
increase with temperature. The former result has been pre- 
viously inferred for x — [Zha et al, 1998; Sinogeikin et al, 
1998], we now show that iron content does not change this 
trend. We also find that iron in the aggregate does not al- 
ter the density discontinuity if its presence does not change 
the volume ratio between phases as it was speculated by 
Yu et al. [2008]. The 0.2-4% density increase estimated by 
a seismic impedance study [Shearer and Flanagan, 1999] is 
consistent with the discontinuity produced in a pyrolite-type 
aggregate, as discussed by Yu et al. [2008]. Overall, iron de- 
creases only slightly magnitudes of contrasts with respect 
to those in the Mg-end member. From our calculations we 
observe that AK < AG and AVp < AVs, which is consis- 
tent with contrasts from PREM [Dziewonski and Anderson, 
1981]. 

4. Conclusions 

For the first time, we have presented parameter-free first- 
principles results of high pressure-temperature aggregate 
elastic properties and sound velocities of Fe-bearing phases 
of olivine and wadsleyite. We used the QHA and a novel 
method of calculating elasticity at high temperatures with- 
out obtaining VDoS of strained configurations [Wu and 
Wentzcovitch, 2011]. Treatment of strain Griineisen param- 
eters via isotropic averages reduced greatly the computa- 
tional cost of the task, which, otherwise, would have been 
prohibitively lengthy. Experiments dealing with simulta- 
neous high pressures and temperatures offer limited data 
and, often relying on ambient conditions, either in temper- 
ature or pressure, to extrapolate to conditions near 410 km 
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(~13.5 GPa and 1500 K). Here, we have shown well de- 
fined changes in elastic and acoustic properties of Fe-bearing 
a— and /?— phases. Contrasts across the a to /3 transition 
are clearly specified near conditions of the 410 km seismic 
discontinuity. Similar sets of results for other coexisting 
phases in the UM and TZ such as pyroxenes, garnets and 
Ca-perovskite are necessary. The consideration of the ef- 
fects of water incorporation in both olivine and wadsleyite 
at relevant conditions of the 410 km discontinuity are also 
critical for this investigation. 
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